GAMLSS_HCP
Calculating centiles based on population-level model.
Kaidi Kang and Simon Vandekar
2022-09-01
knitr::opts_knit$set(root.dir = 'C:/Users/ashua/OneDrive - VUMC/Work/05_OngoingProjects/Lifespan/KaidiReports/GAMLSS_HCP/')
library(gamlss)
library(splines)
library(table1)
library(magrittr)
library(plyr)
library(dplyr)
library(ggplot2)1 Summary
In this project, we would like to implement the GAMLSS model to model the growth curves/charts of the brain using the HCP data.
For this specific report: The deviation score / centiles are calculated based the population-level curves. The variable Family is still included as random effects in the modeling to deal with the correlation between subjects, but when calculating the centiles for each subject, the population-level model (the one without random effects) is used.
2 Data
data = read.csv("HCP_PSYXSEC_DataHarmonized_220829.csv")# some data cleaning
# data$Family[is.na(data$Family)] = paste0("temp", 1:sum(is.na(data$Family)))
# data <- data %>% group_by(Family) %>% slice_sample(n=1)
# nrow(temp) == temp$Family %>% unique %>% length
data$dx_group = data$Diagnosis
data$dx_group[data$Diagnosis != "Normal Control"] = "Schizophrenia Spectrum"
hcp = subset(data, Dataset == "HCP")
# Adding fake family ID to the subject from PSYXSEC
# ind = data$Dataset == "PSYXSEC"
# data$Family[ind] = paste0("f_", 1:sum(ind))
# # subset(data, Dataset == "PSYXSEC")
data$Family %<>% as.factor()
data$group = NA
data$group[data$Dataset == "HCP"] = "Normal - HCP"
data$group[data$dx_group == "Normal Control" & data$Dataset == "PSYXSEC"] = "Normal - PSYXSEC"
data$group[data$dx_group == "Schizophrenia Spectrum"] = "SZ Spectrum - PSYXSEC"Note: The random intercepts/offsets for Family is not included. I noticed that after including the random offset Family, it becomes difficult for the models to converge. For most of the family units in the data, they only have one subjects participating in the studies. Moreover, there are 211 subjects have NA Family ID. This analysis ingored the possible correlation between individuals from the same family.
The random intercepts/offsets for Family are now included in the modeling of the population mean (\(\mu\)) and variance (\(\mu^2 \times \sigma^2\)). All healthy controls from the HCP studies have family ID, however, those from the PSYXSEC study (445 subjects) have missing family ID.
table(hcp$Family) %>% hist(main = "Num of subjects in each family")
# sum(is.na(data$Family[data$Dataset != "HCP"]))table1(~ Age + Site + Race + eTIV + Sex + dx_group | Dataset, data = data)| HCP (N=1952) |
PSYXSEC (N=445) |
Overall (N=2397) |
|
|---|---|---|---|
| Age | |||
| Mean (SD) | 34.5 (19.9) | 29.2 (10.5) | 33.5 (18.7) |
| Median [Min, Max] | 29.0 [5.60, 100] | 25.1 [15.6, 60.5] | 28.0 [5.60, 100] |
| Site | |||
| AMGH | 135 (6.9%) | 0 (0%) | 135 (5.6%) |
| AUCLA | 130 (6.7%) | 0 (0%) | 130 (5.4%) |
| AUMinn | 163 (8.4%) | 0 (0%) | 163 (6.8%) |
| AWashU | 174 (8.9%) | 0 (0%) | 174 (7.3%) |
| DHarvard | 177 (9.1%) | 0 (0%) | 177 (7.4%) |
| DUCLA | 102 (5.2%) | 0 (0%) | 102 (4.3%) |
| DUMinn | 135 (6.9%) | 0 (0%) | 135 (5.6%) |
| DWashU | 111 (5.7%) | 0 (0%) | 111 (4.6%) |
| YUWash | 825 (42.3%) | 0 (0%) | 825 (34.4%) |
| Scan1 | 0 (0%) | 280 (62.9%) | 280 (11.7%) |
| Scan2 | 0 (0%) | 2 (0.4%) | 2 (0.1%) |
| Scan3 | 0 (0%) | 163 (36.6%) | 163 (6.8%) |
| Race | |||
| AA | 272 (13.9%) | 124 (27.9%) | 396 (16.5%) |
| Other | 316 (16.2%) | 23 (5.2%) | 339 (14.1%) |
| White | 1364 (69.9%) | 298 (67.0%) | 1662 (69.3%) |
| eTIV | |||
| Mean (SD) | 1510000 (208000) | 1520000 (193000) | 1510000 (206000) |
| Median [Min, Max] | 1510000 [843000, 2190000] | 1530000 [722000, 2080000] | 1520000 [722000, 2190000] |
| Sex | |||
| F | 1071 (54.9%) | 148 (33.3%) | 1219 (50.9%) |
| M | 881 (45.1%) | 297 (66.7%) | 1178 (49.1%) |
| dx_group | |||
| Normal Control | 1952 (100%) | 211 (47.4%) | 2163 (90.2%) |
| Schizophrenia Spectrum | 0 (0%) | 234 (52.6%) | 234 (9.8%) |
3 GAMLSS
Separate GAMLSS models are fitted to each outcome using the subjects from the HCP study, respectively. Generalized gamma (GG) distribution is specified for each of the outcome.
The final model for the trajectories of all measures are specified as below. We use R.VA as an example:
\[ \mathrm{R.VA} \sim GG(\mu, \sigma, \nu) \]
\[ log(\mu) = \beta_{\mu} + \beta_{\mu, age} \times ns(Age, df = 3) + \beta_{\mu, sex} \times Sex + \beta_{\mu, race} \times Race + \beta_{\mu, eTIV}\times I(neTIV/10000) + \alpha_{Family} \]
\[ log(\sigma) = \beta_{\sigma} + \beta_{\sigma, age} \times ns(Age, df = 3) + \beta_{\sigma, sex} \times Sex + \beta_{\sigma, race} \times Race + \beta_{\sigma, neTIV}\times I(neTIV/10000) + \alpha_{Family} \] \[ \nu = \alpha_{\nu} \]
# control function for gamlss
con_f <- gamlss.control(c.crit = 0.001, n.cyc = 300)
fit_gamlss <- function(y, full.data = data, con = con_f) {
# select the variables needed. The function doesn't allow the existence of missing value in any variables even the variables are not included in the modeling.
# Note: there is a bug in `predict.gamlss` -- when your data frame is called `data`, there will be a problem in predicting...
full.data = full.data[, c(y, "Dataset", "dx_group", "Age", "Sex", "Race", "neTIV", "Family")]
# HCP subjects
hcp = subset(full.data, Dataset == "HCP")
num_obs = nrow(hcp)
hcp = na.omit(hcp)
num_obs_na.rm = nrow(hcp)
cat("Notes:", num_obs - num_obs_na.rm, "observations were removed due to missingness. \n" )
mu.form = as.formula(paste0(y, "~ ns(Age, df = 3) + Sex + Race + I(neTIV/10000) + re(random = ~1|Family)"))
sigma.form = ~ ns(Age, df = 3) + Sex + Race + I(neTIV/10000) + re(random = ~1|Family)
fit = gamlss(mu.form,
sigma.formula = sigma.form,
nu.formula = ~1,
data = hcp,
method = RS(), # or `CG()` or mixed()
family = GG,
control=con)
return(fit)
}# Test
#fit = fit_gamlss("nR.VA")# @param centile_only: if only return centile estimates, no plots. Default = FALSE
plot_gamlss = function(obj, title, full.data = data, centile_only = FALSE){
# 1. DATASET
y = obj$mu.formula[[2]] %>% as.character()
full.data = full.data[, c("Subject", y, "Dataset", "dx_group", "group", "Age", "Sex", "Race", "neTIV", "Family")]
# HCP subjects
hcp = subset(full.data, Dataset == "HCP")
# hcp = model.matrix(obj)
# Schizophrenia Spectrum Subjects
sz = subset(full.data, Dataset != "HCP")
# 2. POPULATION CURVE (averaged over the sample)
if (!centile_only){
age_range = range(hcp$Age)
age_grid = seq(from = age_range[1], to = age_range[2], by = 0.5)
pred_data = hcp[, colnames(hcp) != "Age"]
pred_data = cbind(Age = rep(age_grid, each = nrow(pred_data)), pred_data)
# predict "mu" (without random effects)
# Note: you have to add the original data here....
mean_mu = predict(obj, data = hcp, newdata = pred_data, what = "mu", type = "response")
# predict "sigma"
mean_sigma = predict(obj, data = hcp, newdata = pred_data, what = "sigma", type = "response")
# predict "nu"
mean_nu = predict(obj, data = hcp, newdata = pred_data, what = "nu", type = "response")
## obtain the corresponding 2.5% and 97.5% centiles for mu for each subject
mu_median = gamlss.dist::qGG(0.5, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
mu_LL = gamlss.dist::qGG(0.05, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
mu_UL = gamlss.dist::qGG(0.95, mu = mean_mu, sigma = mean_sigma , nu = mean_nu )
mean_var = mean_sigma^2 * mean_mu^2
pred_data = cbind(mean_mu, mu_median, mu_LL, mu_UL, mean_sigma, mean_var, mean_nu, pred_data)
## marginal means of median & variance at each age value
marg_mean_mu_median = aggregate( mu_median ~ Age*Sex, data = pred_data, mean )
marg_mean_mu_LL = aggregate( mu_LL ~ Age*Sex, data = pred_data, mean )
marg_mean_mu_UL = aggregate( mu_UL ~ Age*Sex, data = pred_data, mean )
marg_mean_sigma = aggregate(mean_sigma ~ Age*Sex, data = pred_data, mean )
marg_mean_var = aggregate(mean_var ~ Age*Sex, data = pred_data, mean)
# 3. PLOT
xlim = c(0, 100)
yrange = range(hcp[,y])
ymin = round_any(yrange[1],10,f=floor)
ymax = round_any(yrange[2],10,f=ceiling)
ylim = c(ymin,ymax)
ybreaks = floor((ylim[2]-ylim[1])/4)
## mu over age
#main = paste0("Population trajectory of the median for ", y, "\n (with 5% and 95% centiles (dotted lines))"),
par(mar = c(4, 5.2, 2.5, 0.5)+0.03)
#,mgp = c(2, 2, 0)) # Set the margin on all sides to 2
colors <- c("pink",
"dodgerblue")
plot(as.formula(paste0(y, "~ Age")), data = hcp, frame.plot = FALSE,
xlab = substitute(paste(bold("Age (years)"))),
xlim = xlim, ylim = ylim,
yaxt = "n",
ylab = expression(bold(paste("Volumes (mm"^3*")"))),
main = title,
col = colors[factor(Sex)],
cex.lab=1.7,
cex.main=2)
#adjustcolor("dodgerblue", alpha.f = 0.35))
axis(side = 2, at=seq(ylim[1],ylim[2],ybreaks))
lines(mu_median ~ Age, data = marg_mean_mu_median[marg_mean_mu_median$Sex=="F",], type = "l", col = "deeppink", lwd = 4)
lines(mu_LL ~ Age, data = marg_mean_mu_LL[marg_mean_mu_median$Sex=="F",], type = "l", lty = 3, col = "deeppink", lwd = 3)
lines(mu_UL ~ Age, data = marg_mean_mu_UL[marg_mean_mu_median$Sex=="F",], type = "l", lty = 3, col = "deeppink", lwd = 3)
lines(mu_median ~ Age, data = marg_mean_mu_median[marg_mean_mu_median$Sex=="M",], type = "l", col = "blue", lwd = 4)
lines(mu_LL ~ Age, data = marg_mean_mu_LL[marg_mean_mu_median$Sex=="M",], type = "l", lty = 3, col = "blue", lwd = 3)
lines(mu_UL ~ Age, data = marg_mean_mu_UL[marg_mean_mu_median$Sex=="M",], type = "l", lty = 3, col = "blue", lwd = 3)
## variance over age
plot(mean_var ~ Age, data = marg_mean_var, type = "l", frame.plot = FALSE,
lwd = 4, col = "dodgerblue",
ylab = "Variance", xlab = "Age",
main = paste0("Normative variance for ", y))
}
# 4. INDIVIDUAL DEVIANCE (i.e., deriving centile conditional on subject characteristics)
# Note: function `centiles()` only works for model fits with one covariate.
# 4.1 predict the centiles for HCP subjects
# ~~sing `predict` function to obtain subject-specific distribution parameter values~~
# ~~Note: the predicted values are based on the values of the random effects.~~
# NEW: the predicted values are NOT based on the random effects
rdata = full.data
## design matrix for HCP subset
ns_obj = ns(hcp$Age, df = 3)
ns_age = predict(ns_obj, rdata$Age)
sex = ifelse(rdata$Sex == "M", 1, 0)
race_other = ifelse(rdata$Race == "Other", 1, 0)
race_white = ifelse(rdata$Race == "White", 1, 0)
neTIV_10000 = rdata$neTIV/10000
mod_mat = cbind(1, ns_age, sex, race_other, race_white, neTIV_10000)
## coef estimates
coef_all = coefAll(obj)
mu_beta_hat = coef_all$mu[1:ncol(mod_mat)]
sigma_beta_hat = coef_all$sigma[1:ncol(mod_mat)]
nu_beta_hat = coef_all$nu
## predicted values
pred_mu = exp(mod_mat %*% mu_beta_hat)
pred_sigma = exp(mod_mat %*% sigma_beta_hat)
pred_nu = rep(nu_beta_hat, nrow(rdata))
rdata = cbind(rdata, pred_mu, pred_sigma, pred_nu)
# pred_mu = predict(obj, data = hcp, what = "mu", type = "response")
# pred_sigma = predict(obj, data = hcp, what = "sigma", type = "response")
# pred_nu = predict(obj, data = hcp, what = "nu", type = "response")
# hcp = cbind(hcp, pred_mu, pred_sigma, pred_nu)
# obtain individual percentiles from subject-specific GG distributions
rdata$subj_centile = gamlss.dist::pGG(rdata[, y], mu = rdata$pred_mu, sigma = rdata$pred_sigma, nu = rdata$pred_nu)
## DENSITY PLOT OF DEVIANCE
if (!centile_only){
plot_centile = ggplot(rdata, aes(x = subj_centile, fill = group, col = group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 30, position="identity") + geom_density(alpha=.2) + labs(title = paste("Centile density plot for", y), x = "Centiles", y = "Density") + facet_wrap(~group, ncol = 1)
centile_med = aggregate( subj_centile ~ group, data = rdata, median)
plot_centile = plot_centile + geom_point(data = centile_med, aes(x = subj_centile, y = 0, col = group), size = 3) + labs(title = y, x = "Centiles", y = "Density")
print(plot_centile)
# # individual deviance (centiles) vs rank
temp_data = subset(rdata, Dataset != "HCP" & dx_group == "Normal Control")
temp_data$rank_n = rank(temp_data$subj_centile)/nrow(temp_data)
plot_centile_v_rank = ggplot(temp_data, aes(x = rank_n, y = subj_centile, col = group)) + geom_point() + geom_abline(slope = 1) + labs(title = paste("rank/sample size vs centiles for", y, "- NC from the PSYXSEC study"), x = "Rank scaled by sample size", y = "Centiles")
print(plot_centile_v_rank)
}
# return centiles
if (centile_only) return(rdata)
}# TEST
# plot_gamlss(fit)4 Model Results
4.1 AV
#fit_R.AV <- fit_gamlss(y = "nR.AV")
#saveRDS(fit_R.AV, file = "fit/fit_R.AV.rds")
fit_R.AV <- readRDS("gamlss_fit_obj/fit_R.AV.rds")
summary(fit_R.AV)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.336486 0.030901 140.334 < 2e-16 ***
## ns(Age, df = 3)1 -0.090390 0.017176 -5.262 1.61e-07 ***
## ns(Age, df = 3)2 -0.395426 0.030133 -13.123 < 2e-16 ***
## ns(Age, df = 3)3 -0.488810 0.023869 -20.479 < 2e-16 ***
## SexM -0.005728 0.007590 -0.755 0.4505
## RaceOther -0.047234 0.011471 -4.118 4.02e-05 ***
## RaceWhite -0.023870 0.009515 -2.509 0.0122 *
## I(neTIV/10000) 0.004947 0.000198 24.989 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.704136 0.169256 -10.068 < 2e-16 ***
## ns(Age, df = 3)1 0.179870 0.084810 2.121 0.03408 *
## ns(Age, df = 3)2 0.094036 0.170651 0.551 0.58168
## ns(Age, df = 3)3 0.210751 0.118277 1.782 0.07496 .
## SexM 0.116232 0.042607 2.728 0.00644 **
## RaceOther -0.067651 0.059252 -1.142 0.25372
## RaceWhite -0.079765 0.048138 -1.657 0.09771 .
## I(neTIV/10000) -0.001974 0.001106 -1.785 0.07448 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.5945 0.3897 4.092 4.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 306.8979
## Residual Deg. of Freedom: 1645.102
## at cycle: 32
##
## Global Deviance: 16991.62
## AIC: 17605.41
## SBC: 19316.86
## ******************************************************************plot_gamlss(fit_R.AV,"Right Anteroventral")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# con_temp <- gamlss.control(c.crit = 0.001, n.cyc = 300, mu.step = 0.1, sigma.step = 0.1, nu.step = 0.1)
# fit_L.AV <- fit_gamlss(y = "nL.AV", con = con_temp)
# saveRDS(fit_L.AV, file = "fit/fit_L.AV.rds")
fit_L.AV <- readRDS("gamlss_fit_obj/fit_L.AV.rds")
summary(fit_L.AV)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3178199 0.0321116 134.463 < 2e-16 ***
## ns(Age, df = 3)1 -0.0987395 0.0170980 -5.775 9.30e-09 ***
## ns(Age, df = 3)2 -0.3806025 0.0308155 -12.351 < 2e-16 ***
## ns(Age, df = 3)3 -0.4929875 0.0231753 -21.272 < 2e-16 ***
## SexM 0.0034004 0.0075974 0.448 0.655
## RaceOther -0.0621468 0.0115845 -5.365 9.35e-08 ***
## RaceWhite -0.0557656 0.0096870 -5.757 1.03e-08 ***
## I(neTIV/10000) 0.0045443 0.0002009 22.623 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.324459 0.162321 -8.159 6.92e-16 ***
## ns(Age, df = 3)1 0.204368 0.083439 2.449 0.0144 *
## ns(Age, df = 3)2 -0.051668 0.151582 -0.341 0.7333
## ns(Age, df = 3)3 0.149183 0.108957 1.369 0.1711
## SexM 0.182792 0.040628 4.499 7.34e-06 ***
## RaceOther -0.079597 0.056698 -1.404 0.1606
## RaceWhite -0.048108 0.046974 -1.024 0.3059
## I(neTIV/10000) -0.004358 0.001049 -4.156 3.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3380 0.3818 8.744 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 412.4777
## Residual Deg. of Freedom: 1539.522
## at cycle: 58
##
## Global Deviance: 16709.26
## AIC: 17534.22
## SBC: 19834.44
## ******************************************************************plot_gamlss(fit_L.AV,"Left Anteroventral")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.2 VA
# fit_R.VA <- fit_gamlss(y = "nR.VA")
# saveRDS(fit_R.VA, file = "fit/fit_R.VA.rds")
fit_R.VA <- readRDS("gamlss_fit_obj/fit_R.VA.rds")
summary(fit_R.VA)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1006038 0.0241937 210.824 < 2e-16 ***
## ns(Age, df = 3)1 -0.1236249 0.0134055 -9.222 < 2e-16 ***
## ns(Age, df = 3)2 -0.3393843 0.0225541 -15.048 < 2e-16 ***
## ns(Age, df = 3)3 -0.4436330 0.0187185 -23.700 < 2e-16 ***
## SexM -0.0104138 0.0060917 -1.710 0.087558 .
## RaceOther -0.0283830 0.0084564 -3.356 0.000809 ***
## RaceWhite 0.0116368 0.0073222 1.589 0.112209
## I(neTIV/10000) 0.0046895 0.0001607 29.187 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.1746573 0.1438867 -15.114 < 2e-16 ***
## ns(Age, df = 3)1 0.2481766 0.0854961 2.903 0.00375 **
## ns(Age, df = 3)2 0.3331504 0.1653229 2.015 0.04406 *
## ns(Age, df = 3)3 0.1901254 0.1165958 1.631 0.10317
## SexM -0.0030207 0.0099475 -0.304 0.76142
## RaceOther -0.0971061 0.0493676 -1.967 0.04936 *
## RaceWhite 0.0109048 0.0169633 0.643 0.52042
## I(neTIV/10000) -0.0013173 0.0007948 -1.657 0.09765 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9922 0.4673 -2.123 0.0339 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 421.4835
## Residual Deg. of Freedom: 1530.517
## at cycle: 33
##
## Global Deviance: 19127.36
## AIC: 19970.32
## SBC: 22320.77
## ******************************************************************fit_R.VA %>% plot_gamlss("Right Ventral Anterior")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.VA <- fit_gamlss(y = "nL.VA")
# saveRDS(fit_L.VA, file = "fit/fit_L.VA.rds")
fit_L.VA <- readRDS("gamlss_fit_obj/fit_L.VA.rds")
summary(fit_L.VA)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0828254 0.0219352 231.721 <2e-16 ***
## ns(Age, df = 3)1 -0.1404496 0.0130618 -10.753 <2e-16 ***
## ns(Age, df = 3)2 -0.4085284 0.0215384 -18.967 <2e-16 ***
## ns(Age, df = 3)3 -0.5232676 0.0209615 -24.963 <2e-16 ***
## SexM 0.0007075 0.0064123 0.110 0.9122
## RaceOther -0.0197369 0.0086920 -2.271 0.0233 *
## RaceWhite 0.0034756 0.0078843 0.441 0.6594
## I(neTIV/10000) 0.0043196 0.0001602 26.964 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3774124 0.1349913 -17.612 < 2e-16 ***
## ns(Age, df = 3)1 0.3204877 0.0857250 3.739 0.000192 ***
## ns(Age, df = 3)2 0.5223207 0.1587604 3.290 0.001026 **
## ns(Age, df = 3)3 0.5748047 0.1119727 5.133 3.23e-07 ***
## SexM 0.0170516 0.0308118 0.553 0.580067
## RaceOther -0.1488782 0.0594528 -2.504 0.012384 *
## RaceWhite -0.1164401 0.0478658 -2.433 0.015109 *
## I(neTIV/10000) -0.0002262 0.0007329 -0.309 0.757665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3881 0.5283 -4.521 6.67e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 494.9939
## Residual Deg. of Freedom: 1457.006
## at cycle: 38
##
## Global Deviance: 18419.8
## AIC: 19409.79
## SBC: 22170.18
## ******************************************************************fit_L.VA %>% plot_gamlss("Left Ventral Anterior")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.3 VL
# fit_R.VL <- fit_gamlss(y = "nR.VL")
# saveRDS(fit_R.VL, file = "fit/fit_R.VL.rds")
fit_R.VL <- readRDS("gamlss_fit_obj/fit_R.VL.rds")
summary(fit_R.VL)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.431e+00 1.631e-02 394.423 < 2e-16 ***
## ns(Age, df = 3)1 -1.062e-01 1.130e-02 -9.397 < 2e-16 ***
## ns(Age, df = 3)2 -3.341e-01 1.695e-02 -19.709 < 2e-16 ***
## ns(Age, df = 3)3 -4.270e-01 1.521e-02 -28.063 < 2e-16 ***
## SexM -3.651e-03 1.171e-03 -3.119 0.00185 **
## RaceOther -1.442e-02 6.148e-03 -2.345 0.01919 *
## RaceWhite 8.041e-03 5.643e-03 1.425 0.15439
## I(neTIV/10000) 3.936e-03 9.319e-05 42.235 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.1289540 0.1243124 -17.126 < 2e-16 ***
## ns(Age, df = 3)1 0.1787333 0.1855895 0.963 0.335693
## ns(Age, df = 3)2 0.5052426 0.2433606 2.076 0.038074 *
## ns(Age, df = 3)3 0.4293012 0.1866583 2.300 0.021605 *
## SexM 0.0109421 0.0117163 0.934 0.350513
## RaceOther -0.1572408 0.0682413 -2.304 0.021364 *
## RaceWhite -0.0956725 0.0350703 -2.728 0.006455 **
## I(neTIV/10000) -0.0030741 0.0008991 -3.419 0.000648 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.29756 0.06088 -4.888 1.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 611.083
## Residual Deg. of Freedom: 1340.917
## at cycle: 31
##
## Global Deviance: 22859.79
## AIC: 24081.95
## SBC: 27489.73
## ******************************************************************fit_R.VL %>% plot_gamlss("Right Ventrolateral")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.VL <- fit_gamlss(y = "nL.VL")
# saveRDS(fit_L.VL, file = "fit/fit_L.VL.rds")
fit_L.VL <- readRDS("gamlss_fit_obj/fit_L.VL.rds")
summary(fit_L.VL)## Warning in summary.gamlss(fit_L.VL): summary: vcov has failed, option qr is used instead## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,
## family = GG, data = hcp, method = RS(), control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.4164555 0.0180096 356.280 <2e-16 ***
## ns(Age, df = 3)1 -0.1255590 0.0101464 -12.375 <2e-16 ***
## ns(Age, df = 3)2 -0.4018762 0.0166799 -24.093 <2e-16 ***
## ns(Age, df = 3)3 -0.4768878 0.0145035 -32.881 <2e-16 ***
## SexM -0.0070763 0.0045983 -1.539 0.1240
## RaceOther -0.0167174 0.0066761 -2.504 0.0124 *
## RaceWhite -0.0038386 0.0057564 -0.667 0.5050
## I(neTIV/10000) 0.0039653 0.0001202 32.993 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.168233 0.163836 -13.234 < 2e-16 ***
## ns(Age, df = 3)1 0.341323 0.085429 3.995 6.79e-05 ***
## ns(Age, df = 3)2 0.545283 0.161556 3.375 0.000758 ***
## ns(Age, df = 3)3 0.386741 0.115105 3.360 0.000800 ***
## SexM 0.039165 0.041883 0.935 0.349893
## RaceOther -0.132506 0.059284 -2.235 0.025566 *
## RaceWhite -0.087382 0.048322 -1.808 0.070767 .
## I(neTIV/10000) -0.003180 0.001075 -2.957 0.003153 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5194 0.3969 -1.308 0.191
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 540.4578
## Residual Deg. of Freedom: 1411.542
## at cycle: 30
##
## Global Deviance: 22688.86
## AIC: 23769.77
## SBC: 26783.7
## ******************************************************************fit_L.VL %>% plot_gamlss("Left Ventrolateral")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.4 VPL
# con_temp <- gamlss.control(c.crit = 0.001, n.cyc = 300, mu.step = 0.01, sigma.step = 0.01, nu.step = 0.01)
# fit_R.VPL <- fit_gamlss(y = "nR.VPL", con = con_temp)
# saveRDS(fit_R.VPL, file = "fit/fit_R.VPL.rds")
fit_R.VPL <- readRDS("gamlss_fit_obj/fit_R.VPL.rds")
summary(fit_R.VPL)## Warning in summary.gamlss(fit_R.VPL): summary: vcov has failed, option qr is used instead## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,
## family = GG, data = hcp, method = RS(), control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1381379 0.0215518 238.409 < 2e-16 ***
## ns(Age, df = 3)1 0.0024010 0.0116605 0.206 0.8369
## ns(Age, df = 3)2 -0.0807364 0.0198881 -4.060 5.11e-05 ***
## ns(Age, df = 3)3 -0.1880015 0.0163009 -11.533 < 2e-16 ***
## SexM 0.0340769 0.0054827 6.215 6.25e-10 ***
## RaceOther 0.0477738 0.0080571 5.929 3.59e-09 ***
## RaceWhite 0.0150638 0.0068230 2.208 0.0274 *
## I(neTIV/10000) 0.0032241 0.0001426 22.616 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.234192 0.163893 -13.632 <2e-16 ***
## ns(Age, df = 3)1 0.188653 0.085462 2.207 0.0274 *
## ns(Age, df = 3)2 0.555743 0.161604 3.439 0.0006 ***
## ns(Age, df = 3)3 0.269718 0.115153 2.342 0.0193 *
## SexM 0.040075 0.041898 0.956 0.3390
## RaceOther -0.094750 0.059306 -1.598 0.1103
## RaceWhite -0.071702 0.048342 -1.483 0.1382
## I(neTIV/10000) -0.001886 0.001076 -1.754 0.0797 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01304 0.34109 -0.038 0.969
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 500.9787
## Residual Deg. of Freedom: 1451.021
## at cycle: 36
##
## Global Deviance: 18503.05
## AIC: 19505.01
## SBC: 22298.77
## ******************************************************************fit_R.VPL %>% plot_gamlss("Right Ventral Posterolateral")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
#
# fit_L.VPL <- fit_gamlss(y = "nL.VPL")
# saveRDS(fit_L.VPL, file = "fit/fit_L.VPL.rds")
fit_L.VPL <- readRDS("gamlss_fit_obj/fit_L.VPL.rds")
summary(fit_L.VPL)## Warning in summary.gamlss(fit_L.VPL): summary: vcov has failed, option qr is used instead## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,
## family = GG, data = hcp, method = RS(), control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.175952 0.022239 232.740 < 2e-16 ***
## ns(Age, df = 3)1 0.002724 0.012273 0.222 0.8244
## ns(Age, df = 3)2 -0.112207 0.020710 -5.418 6.78e-08 ***
## ns(Age, df = 3)3 -0.220563 0.017225 -12.805 < 2e-16 ***
## SexM 0.026222 0.005690 4.608 4.33e-06 ***
## RaceOther 0.032999 0.008381 3.937 8.53e-05 ***
## RaceWhite 0.016827 0.007156 2.351 0.0188 *
## I(neTIV/10000) 0.003234 0.000148 21.855 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.222516 0.163864 -13.563 < 2e-16 ***
## ns(Age, df = 3)1 0.199709 0.085445 2.337 0.019552 *
## ns(Age, df = 3)2 0.540559 0.161578 3.346 0.000841 ***
## ns(Age, df = 3)3 0.298066 0.115129 2.589 0.009716 **
## SexM 0.051929 0.041890 1.240 0.215294
## RaceOther -0.117022 0.059294 -1.974 0.048608 *
## RaceWhite -0.081936 0.048332 -1.695 0.090221 .
## I(neTIV/10000) -0.001655 0.001076 -1.539 0.124115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3206 0.3367 -0.952 0.341
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 414.7198
## Residual Deg. of Freedom: 1537.28
## at cycle: 36
##
## Global Deviance: 18722.16
## AIC: 19551.6
## SBC: 21864.33
## ******************************************************************fit_L.VPL %>% plot_gamlss("Left Ventral Posterolateral")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.5 Pul
# fit_R.Pul <- fit_gamlss(y = "nR.Pul")
# saveRDS(fit_R.Pul, file = "fit/fit_R.Pul.rds")
fit_R.Pul <- readRDS("gamlss_fit_obj/fit_R.Pul.rds")
summary(fit_R.Pul)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9999426 0.0198245 353.095 < 2e-16 ***
## ns(Age, df = 3)1 -0.1540033 0.0113824 -13.530 < 2e-16 ***
## ns(Age, df = 3)2 -0.5157236 0.0184744 -27.916 < 2e-16 ***
## ns(Age, df = 3)3 -0.5236371 0.0165736 -31.595 < 2e-16 ***
## SexM 0.0357509 0.0051067 7.001 3.9e-12 ***
## RaceOther 0.0695703 0.0076612 9.081 < 2e-16 ***
## RaceWhite 0.0634289 0.0064723 9.800 < 2e-16 ***
## I(neTIV/10000) 0.0018047 0.0001318 13.689 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.7560376 0.1433865 -19.221 < 2e-16 ***
## ns(Age, df = 3)1 0.2710613 0.0835099 3.246 0.001198 **
## ns(Age, df = 3)2 0.6899750 0.1600834 4.310 1.74e-05 ***
## ns(Age, df = 3)3 0.4070112 0.1094739 3.718 0.000209 ***
## SexM 0.0073490 0.0250248 0.294 0.769055
## RaceOther -0.0989740 0.0592871 -1.669 0.095256 .
## RaceWhite -0.1359840 0.0475934 -2.857 0.004336 **
## I(neTIV/10000) 0.0012873 0.0008579 1.500 0.133715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3747 0.5348 6.31 3.72e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 519.4501
## Residual Deg. of Freedom: 1432.55
## at cycle: 31
##
## Global Deviance: 24211.35
## AIC: 25250.25
## SBC: 28147.02
## ******************************************************************fit_R.Pul %>% plot_gamlss("Right Pulvinar")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.Pul <- fit_gamlss(y = "nL.Pul")
# saveRDS(fit_L.Pul, file = "fit/fit_L.Pul.rds")
fit_L.Pul <- readRDS("gamlss_fit_obj/fit_L.Pul.rds")
summary(fit_L.Pul)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9976200 0.0203387 344.055 < 2e-16 ***
## ns(Age, df = 3)1 -0.1386545 0.0119464 -11.606 < 2e-16 ***
## ns(Age, df = 3)2 -0.4880280 0.0197342 -24.730 < 2e-16 ***
## ns(Age, df = 3)3 -0.4868830 0.0191067 -25.482 < 2e-16 ***
## SexM 0.0330419 0.0052252 6.324 3.34e-10 ***
## RaceOther 0.0749004 0.0078438 9.549 < 2e-16 ***
## RaceWhite 0.0773936 0.0065637 11.791 < 2e-16 ***
## I(neTIV/10000) 0.0020230 0.0001347 15.018 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.403251 0.162671 -14.774 < 2e-16 ***
## ns(Age, df = 3)1 0.238463 0.084560 2.820 0.00486 **
## ns(Age, df = 3)2 0.809308 0.159989 5.059 4.73e-07 ***
## ns(Age, df = 3)3 0.684229 0.113063 6.052 1.80e-09 ***
## SexM 0.041448 0.042030 0.986 0.32422
## RaceOther -0.066980 0.060658 -1.104 0.26967
## RaceWhite -0.110740 0.048557 -2.281 0.02271 *
## I(neTIV/10000) -0.001283 0.001082 -1.186 0.23571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7255 0.5318 5.126 3.34e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 412.4122
## Residual Deg. of Freedom: 1539.588
## at cycle: 17
##
## Global Deviance: 24444.56
## AIC: 25269.38
## SBC: 27569.24
## ******************************************************************fit_L.Pul %>% plot_gamlss("Left Pulvinar")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.6 LGN
# fit_R.LGN <- fit_gamlss(y = "nR.LGN")
# saveRDS(fit_R.LGN, file = "fit/fit_R.LGN.rds")
fit_R.LGN <- readRDS("gamlss_fit_obj/fit_R.LGN.rds")
summary(fit_R.LGN)## Warning in summary.gamlss(fit_R.LGN): summary: vcov has failed, option qr is used instead## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,
## family = GG, data = hcp, method = RS(), control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.2011170 0.0257284 163.287 < 2e-16 ***
## ns(Age, df = 3)1 -0.0555857 0.0146651 -3.790 0.000155 ***
## ns(Age, df = 3)2 -0.3102835 0.0249621 -12.430 < 2e-16 ***
## ns(Age, df = 3)3 -0.5386605 0.0232848 -23.134 < 2e-16 ***
## SexM 0.0178171 0.0066408 2.683 0.007359 **
## RaceOther 0.0015038 0.0094978 0.158 0.874211
## RaceWhite -0.0371695 0.0081972 -4.534 6.13e-06 ***
## I(neTIV/10000) 0.0025200 0.0001721 14.643 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.3151720 0.1638457 -14.130 < 2e-16 ***
## ns(Age, df = 3)1 0.1182908 0.0854323 1.385 0.16637
## ns(Age, df = 3)2 0.6428284 0.1615588 3.979 7.25e-05 ***
## ns(Age, df = 3)3 0.5774035 0.1150985 5.017 5.87e-07 ***
## SexM -0.0152082 0.0418851 -0.363 0.71659
## RaceOther -0.1596555 0.0592881 -2.693 0.00716 **
## RaceWhite -0.0909612 0.0483254 -1.882 0.05999 .
## I(neTIV/10000) 0.0003481 0.0010754 0.324 0.74619
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3456 0.2830 1.221 0.222
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 431.9813
## Residual Deg. of Freedom: 1520.019
## at cycle: 32
##
## Global Deviance: 14754.76
## AIC: 15618.72
## SBC: 18027.71
## ******************************************************************fit_R.LGN %>% plot_gamlss("Right Lateral Geniculate")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.LGN <- fit_gamlss(y = "nL.LGN")
# saveRDS(fit_L.LGN, file = "fit/fit_L.LGN.rds")
fit_L.LGN <- readRDS("gamlss_fit_obj/fit_L.LGN.rds")
summary(fit_L.LGN)## Warning in summary.gamlss(fit_L.LGN): summary: vcov has failed, option qr is used instead## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form, nu.formula = ~1,
## family = GG, data = hcp, method = RS(), control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.1886188 0.0268914 155.761 < 2e-16 ***
## ns(Age, df = 3)1 -0.0706637 0.0154386 -4.577 5.01e-06 ***
## ns(Age, df = 3)2 -0.3468106 0.0254375 -13.634 < 2e-16 ***
## ns(Age, df = 3)3 -0.5853796 0.0222919 -26.260 < 2e-16 ***
## SexM 0.0369125 0.0069437 5.316 1.18e-07 ***
## RaceOther 0.0126113 0.0101088 1.248 0.2123
## RaceWhite -0.0138185 0.0083933 -1.646 0.0998 .
## I(neTIV/10000) 0.0025617 0.0001792 14.292 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.238e+00 1.638e-01 -13.658 < 2e-16 ***
## ns(Age, df = 3)1 3.200e-01 8.542e-02 3.747 0.000186 ***
## ns(Age, df = 3)2 4.770e-01 1.615e-01 2.953 0.003196 **
## ns(Age, df = 3)3 3.497e-01 1.151e-01 3.038 0.002417 **
## SexM -1.818e-02 4.188e-02 -0.434 0.664297
## RaceOther -3.015e-02 5.928e-02 -0.509 0.611107
## RaceWhite -5.679e-02 4.832e-02 -1.175 0.240041
## I(neTIV/10000) -2.055e-05 1.075e-03 -0.019 0.984758
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4187 0.2778 -1.507 0.132
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 363.5109
## Residual Deg. of Freedom: 1588.489
## at cycle: 20
##
## Global Deviance: 14953.17
## AIC: 15680.19
## SBC: 17707.35
## ******************************************************************fit_L.LGN %>% plot_gamlss("Left Lateral Geniculate")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.7 MGN
# fit_R.MGN <- fit_gamlss(y = "nR.MGN")
# saveRDS(fit_R.MGN, file = "fit/fit_R.MGN.rds")
fit_R.MGN <- readRDS("gamlss_fit_obj/fit_R.MGN.rds")
summary(fit_R.MGN)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8486749 0.0281790 136.580 < 2e-16 ***
## ns(Age, df = 3)1 -0.1012390 0.0157224 -6.439 1.56e-10 ***
## ns(Age, df = 3)2 -0.5185034 0.0271597 -19.091 < 2e-16 ***
## ns(Age, df = 3)3 -0.5793734 0.0232749 -24.893 < 2e-16 ***
## SexM 0.0300637 0.0072555 4.144 3.59e-05 ***
## RaceOther 0.0527740 0.0106932 4.935 8.79e-07 ***
## RaceWhite 0.0368115 0.0088847 4.143 3.59e-05 ***
## I(neTIV/10000) 0.0020486 0.0001869 10.961 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.2799009 0.1522172 -14.978 < 2e-16 ***
## ns(Age, df = 3)1 0.1152035 0.0816016 1.412 0.158199
## ns(Age, df = 3)2 0.4628272 0.1615942 2.864 0.004233 **
## ns(Age, df = 3)3 0.3692912 0.1083389 3.409 0.000668 ***
## SexM -0.0164405 0.0357480 -0.460 0.645648
## RaceOther -0.0524711 0.0581812 -0.902 0.367262
## RaceWhite -0.0822500 0.0474837 -1.732 0.083425 .
## I(neTIV/10000) 0.0006995 0.0009383 0.745 0.456101
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9616 0.4276 4.587 4.82e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 262.0145
## Residual Deg. of Freedom: 1689.986
## at cycle: 26
##
## Global Deviance: 13187.95
## AIC: 13711.98
## SBC: 15173.13
## ******************************************************************fit_R.MGN %>% plot_gamlss("Right Medial Geniculate")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.MGN <- fit_gamlss(y = "nL.MGN")
# saveRDS(fit_L.MGN, file = "fit/fit_L.MGN.rds")
fit_L.MGN <- readRDS("gamlss_fit_obj/fit_L.MGN.rds")
summary(fit_L.MGN)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8313748 0.0262356 146.037 < 2e-16 ***
## ns(Age, df = 3)1 -0.0993767 0.0146293 -6.793 1.55e-11 ***
## ns(Age, df = 3)2 -0.4568395 0.0248393 -18.392 < 2e-16 ***
## ns(Age, df = 3)3 -0.5206643 0.0202392 -25.726 < 2e-16 ***
## SexM 0.0251153 0.0067533 3.719 0.000207 ***
## RaceOther 0.0425273 0.0099682 4.266 2.10e-05 ***
## RaceWhite 0.0196793 0.0080555 2.443 0.014675 *
## I(neTIV/10000) 0.0020520 0.0001746 11.756 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.4228414 0.1619408 -14.961 < 2e-16 ***
## ns(Age, df = 3)1 0.2084008 0.0837955 2.487 0.01298 *
## ns(Age, df = 3)2 0.4420514 0.1620983 2.727 0.00646 **
## ns(Age, df = 3)3 0.2158756 0.1135599 1.901 0.05748 .
## SexM -0.0083584 0.0482355 -0.173 0.86245
## RaceOther 0.0185452 0.0743993 0.249 0.80319
## RaceWhite -0.0457813 0.0558899 -0.819 0.41283
## I(neTIV/10000) 0.0008301 0.0011769 0.705 0.48071
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6791 0.4233 3.966 7.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 358.2205
## Residual Deg. of Freedom: 1593.78
## at cycle: 13
##
## Global Deviance: 12902.2
## AIC: 13618.65
## SBC: 15616.3
## ******************************************************************fit_L.MGN %>% plot_gamlss("Left Medial Geniculate")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.8 CM
# fit_R.CM <- fit_gamlss(y = "nR.CM")
# saveRDS(fit_R.CM, file = "fit/fit_R.CM.rds")
fit_R.CM <- readRDS("gamlss_fit_obj/fit_R.CM.rds")
summary(fit_R.CM)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4991495 0.0311324 144.517 < 2e-16 ***
## ns(Age, df = 3)1 -0.0859788 0.0179499 -4.790 1.82e-06 ***
## ns(Age, df = 3)2 -0.4876218 0.0299364 -16.289 < 2e-16 ***
## ns(Age, df = 3)3 -0.4529192 0.0255408 -17.733 < 2e-16 ***
## SexM 0.0391498 0.0081736 4.790 1.82e-06 ***
## RaceOther 0.0890777 0.0116438 7.650 3.37e-14 ***
## RaceWhite 0.0657958 0.0096864 6.793 1.53e-11 ***
## I(neTIV/10000) 0.0015034 0.0002095 7.177 1.07e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.531210 0.166052 -15.243 < 2e-16 ***
## ns(Age, df = 3)1 0.258450 0.082778 3.122 0.00183 **
## ns(Age, df = 3)2 0.492311 0.160946 3.059 0.00226 **
## ns(Age, df = 3)3 0.249255 0.108104 2.306 0.02125 *
## SexM -0.058628 0.042169 -1.390 0.16462
## RaceOther -0.040915 0.058948 -0.694 0.48772
## RaceWhite -0.046447 0.047547 -0.977 0.32878
## I(neTIV/10000) 0.002861 0.001079 2.651 0.00811 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.2473 0.3835 5.86 5.57e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 288.2003
## Residual Deg. of Freedom: 1663.8
## at cycle: 24
##
## Global Deviance: 15969.16
## AIC: 16545.56
## SBC: 18152.74
## ******************************************************************fit_R.CM %>% plot_gamlss("Right Centromedian")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.CM <- fit_gamlss(y = "nL.CM")
# saveRDS(fit_L.CM, file = "fit/fit_L.CM.rds")
fit_L.CM <- readRDS("gamlss_fit_obj/fit_L.CM.rds")
summary(fit_L.CM)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5114604 0.0285677 157.922 < 2e-16 ***
## ns(Age, df = 3)1 -0.0729940 0.0165206 -4.418 1.06e-05 ***
## ns(Age, df = 3)2 -0.4496853 0.0268738 -16.733 < 2e-16 ***
## ns(Age, df = 3)3 -0.4193886 0.0228457 -18.357 < 2e-16 ***
## SexM 0.0362262 0.0074598 4.856 1.32e-06 ***
## RaceOther 0.0701914 0.0114136 6.150 9.84e-10 ***
## RaceWhite 0.0669973 0.0092582 7.237 7.20e-13 ***
## I(neTIV/10000) 0.0015892 0.0001913 8.308 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.555025 0.162582 -15.715 < 2e-16 ***
## ns(Age, df = 3)1 0.298060 0.082560 3.610 0.000316 ***
## ns(Age, df = 3)2 0.554606 0.162581 3.411 0.000663 ***
## ns(Age, df = 3)3 0.258942 0.111402 2.324 0.020233 *
## SexM -0.049854 0.040910 -1.219 0.223171
## RaceOther 0.001783 0.011640 0.153 0.878302
## RaceWhite -0.118319 0.036501 -3.242 0.001214 **
## I(neTIV/10000) 0.002594 0.001074 2.415 0.015852 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9395 0.3757 5.163 2.74e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 394.9273
## Residual Deg. of Freedom: 1557.073
## at cycle: 26
##
## Global Deviance: 15821.92
## AIC: 16611.77
## SBC: 18814.13
## ******************************************************************fit_L.CM %>% plot_gamlss("Left Centromedian")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
4.9 MD
# fit_R.MD <- fit_gamlss(y = "nR.MD")
# saveRDS(fit_R.MD, file = "fit/fit_R.MD.rds")
fit_R.MD <- readRDS("gamlss_fit_obj/fit_R.MD.rds")
summary(fit_R.MD)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1509257 0.0217291 283.074 < 2e-16 ***
## ns(Age, df = 3)1 -0.0534377 0.0130422 -4.097 4.40e-05 ***
## ns(Age, df = 3)2 -0.3788348 0.0218202 -17.362 < 2e-16 ***
## ns(Age, df = 3)3 -0.3209589 0.0195392 -16.426 < 2e-16 ***
## SexM -0.0055777 0.0059875 -0.932 0.352
## RaceOther 0.0625791 0.0083124 7.528 8.84e-14 ***
## RaceWhite 0.0500609 0.0068067 7.355 3.14e-13 ***
## I(neTIV/10000) 0.0015626 0.0001506 10.377 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.1913203 0.1482312 -21.529 < 2e-16 ***
## ns(Age, df = 3)1 0.1317723 0.0744351 1.770 0.07688 .
## ns(Age, df = 3)2 0.8101773 0.1514237 5.350 1.01e-07 ***
## ns(Age, df = 3)3 0.4673210 0.0980506 4.766 2.06e-06 ***
## SexM -0.0982921 0.0369703 -2.659 0.00793 **
## RaceOther 0.0054464 0.0291271 0.187 0.85169
## RaceWhite -0.0229580 0.0359510 -0.639 0.52319
## I(neTIV/10000) 0.0042454 0.0009619 4.414 1.09e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.5742 0.8903 10.75 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 457.7676
## Residual Deg. of Freedom: 1494.232
## at cycle: 76
##
## Global Deviance: 21552.63
## AIC: 22468.16
## SBC: 25020.95
## ******************************************************************fit_R.MD %>% plot_gamlss("Right Mediodorsal")## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "mu", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
# fit_L.MD <- fit_gamlss(y = "nL.MD")
# saveRDS(fit_L.MD, file = "fit/fit_L.MD.rds")
fit_L.MD <- readRDS("gamlss_fit_obj/fit_L.MD.rds")
summary(fit_L.MD)## ******************************************************************
## Family: c("GG", "generalised Gamma Lopatatsidis-Green")
##
## Call: gamlss(formula = mu.form, sigma.formula = sigma.form,
## nu.formula = ~1, family = GG, data = hcp, method = RS(),
## control = con)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1434077 0.0277858 221.099 < 2e-16 ***
## ns(Age, df = 3)1 -0.0438666 0.0159427 -2.752 0.006 **
## ns(Age, df = 3)2 -0.4067019 0.0267378 -15.211 < 2e-16 ***
## ns(Age, df = 3)3 -0.3461348 0.0227581 -15.209 < 2e-16 ***
## SexM 0.0024296 0.0071370 0.340 0.734
## RaceOther 0.0508062 0.0106681 4.762 2.09e-06 ***
## RaceWhite 0.0500736 0.0084953 5.894 4.60e-09 ***
## I(neTIV/10000) 0.0014357 0.0001861 7.714 2.17e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.558733 0.159370 -16.055 < 2e-16 ***
## ns(Age, df = 3)1 0.164944 0.083612 1.973 0.04870 *
## ns(Age, df = 3)2 0.488547 0.162819 3.001 0.00274 **
## ns(Age, df = 3)3 0.308737 0.110398 2.797 0.00523 **
## SexM -0.046958 0.040568 -1.158 0.24724
## RaceOther 0.040273 0.058324 0.691 0.48997
## RaceWhite -0.056568 0.047855 -1.182 0.23736
## I(neTIV/10000) 0.002360 0.001054 2.239 0.02532 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.942 0.432 6.811 1.38e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1952
## Degrees of Freedom for the fit: 393.6556
## Residual Deg. of Freedom: 1558.344
## at cycle: 34
##
## Global Deviance: 21899.1
## AIC: 22686.41
## SBC: 24881.68
## ******************************************************************fit_L.MD %>% plot_gamlss("Left Mediodorsal")## Warning in data.frame(..., check.names = FALSE): row names were found from a short variable and have been discarded
## Warning in data.frame(..., check.names = FALSE): There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction## Warning in predict.gamlss(obj, data = hcp, newdata = pred_data, what = "sigma", : There is a discrepancy between the original and the re-fit
## used to achieve 'safe' predictions
## ## new prediction
5 Centile Mahalanobis Distance (CMD)
To create a composite score to indicate an individual’s overall deviation from the normative curve, CMD is calculated for each subject. CMD for the \(i\)-th subject is formalized as \[ D_i(x) = \sqrt{(x - \mu)^T S^{-1} (x-\mu) } \]
where \(x\) is the obtained centiles for each subject, \(\mu\) and \(S\) are the means and covariance of centiles among the normal controls (HCP subjects).
# adding predicted centiles with respect to each outcome to the data
outcome_var = paste0(rep(c("nR.", "nL."), times = 9), rep(c("AV", "VA", "VL", "VPL", "Pul", "LGN", "MGN", "CM", "MD"), each = 2))
output_data = data[, c("Subject", "Family", "Diagnosis", "dx_group", "Dataset", "group", "Age", "Sex", "Race", "Site", "neTIV", outcome_var)]
# obtaining the centile for each outcome
for (k in outcome_var){
outcome = substring(k, 2)
fit = get(paste0("fit_", outcome))
centiles = plot_gamlss(fit, centile_only = TRUE)
# subject centile for this outcome
output_data = output_data[order(output_data$Subject), ]
centiles = centiles[order(centiles$Subject), ]
output_data = merge(output_data, centiles[, c("Subject", "subj_centile")], by = "Subject")
colnames(output_data)[colnames(output_data) == "subj_centile"] = paste0("centile_", outcome)
}
# CMD
ind_row = output_data$group == "Normal - HCP"
ind_col = substring(colnames(output_data), 1, 8) == "centile_"
S = cov(output_data[ind_row, ind_col])
mu = colMeans(output_data[ind_row, ind_col])
center_X = apply(output_data[, ind_col], 1, function(y) {y - mu}) %>% t()
CMD = apply(center_X, 1, function(y) sqrt(c(y) %*% solve(S) %*% c(y)))
output_data = cbind(output_data, CMD)write.csv(output_data, "data_centiles_main.csv")plot_CMD = ggplot(output_data, aes(x = CMD, fill = group, col = group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 50, position="identity") + geom_density(alpha=.2) + labs(x = "Centiles", y = "Density")
# + facet_wrap(~group, ncol = 1)
CMD_med = aggregate( CMD ~ group, data = output_data, median)
plot_CMD = plot_CMD + geom_point(data = CMD_med, aes(x = CMD, y = 0, col = group), size = 3) + labs(title = "Density of CMD by groups", x = "Centiles", y = "Density")
print(plot_CMD)
plot_CMD = ggplot(output_data, aes(x = CMD, fill = dx_group, col = dx_group)) + geom_histogram(aes(y=..density..), alpha = 0.2, bins = 50, position="identity") + geom_density(alpha=.2) + labs(x = "Centiles", y = "Density")
# + facet_wrap(~group, ncol = 1)
CMD_med = aggregate( CMD ~ dx_group, data = output_data, median)
plot_CMD = plot_CMD + geom_point(data = CMD_med, aes(x = CMD, y = 0, col = dx_group), size = 3) + labs(title = "Density of CMD by diagnosis", x = "Centiles", y = "Density")
print(plot_CMD)